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Abstract. The influence of non magnetic impurities in the 2d XY model is investigated through Monte 
Carlo (MC) simulations. The general picture of the transition is fully understood from the Harris criterion 
which predicts that the universality class is unchanged, and the Berezinskii-Kosterlitz-Thouless description 
of the topological transition remains valid. We nevertheless address here the question about the influence 
of dilution on the quasi-long-range order at low temperatures. In particular, we study the asymptotic of the 
pair correlation function and report the MC estimates for the critical exponent rj at different dilutions. In 
the weak dilution region, our MC calculations are further supported by simple spin-wave-like calculations. 

Key words. XY model - topological transition - random systems. 

PACS. 05.50.+q Lattice theory and statistics (Ising, Potts, etc.) - 64.60.Fr Equilibrium properties near 
critical points, critical exponents - 75.10.Hk Classical spin models 

1 Introduction sumes that the disorientation 9r — 6'r-f/x between neigh- 
bouring spins {fi denotes the unit lattice spacing) remains 

From general scaling arguments, Harris criterion predicts small, captures the essential features of the behaviour of 

that Gaussian disorder (coupled to the energy density) is the system, leading to an algebraic decay of the correlation 

an irrelevant perturbation at the fixed point of any pure function [7] 
system when the exponent a which governs the singular 

behaviour of the specific heat of the pure model in the (""(r) •cr(r + R)) ~ |R|^'''"^\ (1) 
vicinity of the critical point is negative [1]. The critical 

behaviour of the random system is thus governed by the (where cr(r) are two-component unit vectors), and, as a 

pure fixed point and the universality class of the model consequence, to the absence of magnetization and to an 

remains unchanged. infinite susceptibility. The correlation function exponent. 

From this respect, the 2c? XY model is probably the 77(T), continuously increases with temperature. The model 

one where the influence of quenched disorder is essentially is said to have a low temperature phase with quasi-long- 

trivial, and nothing special is expected, since it displays a range order. At high temperatures on the other hand, a 

critical behaviour described by essential singularities. On high temperature series expansion leads to a more con- 

the other hand, the transition of the 2d XY model is not ventional exponential decay of the correlation function. It 

conventional and has very interesting features which make is clear from these two extreme behaviours that something 

the question of the influence of disorder worth studying, must happen in between and that the model undergoes a 

Indeed, the transition, first described by Berezinskii [2] transition in the intermediate regime. The scenario pro- 

and Kosterlitz and Thouless [3,4], is governed by the con- posed by Kosterlitz and Thouless (KT) is based on the 

densation of topological defects. The symmetry of the existence of vortices, which are localized defects where 

model being continuous, there is no spontaneous order in the field 9^ (the angle between spin (T(r) and some ar- 

the system (in 2d) at any non zero temperature accord- bitrary reference direction) may become singular at some 

ing to Mermin Wagner theorem [5,6]. Nevertheless at very points. The energy carried by such defects increases log- 

low temperatures, the spin-wave approximation which as- arithmically with the system size, and thus they are con- 
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strained to be associated in pairs which nevertheless am- 
plify the disordering of the system, leading to an effec- 
tive increase of the temperature. The vortices appear in 
increasing number in the system when the temperature 
increases, and the transition temperature is reached when 
the pairs break, leaving the system totally disordered. In 
the high temperature phase, the approach to the critical 
point Tc, (also referred to as KT temperature, Tkt), is 
described by essential singularities, for instance, the cor- 
relation length and the susceptibility behave as ^(i) ~ 
exp(-6t-'^), ~ ^^"^ t=\T- TktI, with 77 = i and 
a= \ (for reviews, sec e.g. [8,9,10,11,12]). 

From the hyperscaling relation a = 2 — vd, and from 
the fact that the correlation length has an essential singu- 
larity at the critical point, one deduces a non typical value 
for the exponent of the specific heat, a = —00, which im- 
plies, according to Harris criterion, that the adjunction of 
thermal randomness (i.e. coupled to the energy density) 
does not modify the general scenario described above. All 
along the transition line, the KT universality class might 
be recovered (that is an exponent r7(TKT) = \ for the 
correlation function decay). Let us specify the case of site 
dilution for example. The original Hamiltonian of the pure 
2d XY model, 

r n 

where cr(r) • cr(r -I- /i) denotes the scalar product and the 
spins cr(r) are located at the lattice sites r of a square lat- 
tice A, is modified by the introduction of a set of occupa- 
tion variables Cr which take the values (with probability 
1 — p) or 1 (with probability p) depending on the fact that 
site r is empty or occupied by a spin, 

P[Cr]=n^(Cr) 
r 

= n[P^(Cr-l) + (l-p)5(Cr)]. (3) 
r 

Thus the Hamiltonian of the diluted 2d XY model reads: 
H = -jY^Y^CrCr+^cr{Y)-a{Y + n). (4) 

r M 

Altogether, one expects a phase diagram starting from the 
pure system critical temperature fc^TxT/^ — 0.893 [13] at 
p = 1, and decreasing up to a transition at zero temper- 
ature at the site percolation threshold of the 2d square 
lattice, pc — 0.59, since no transition at all takes place 
in the system when there is no more percolating cluster 
of spins. The question which remains interesting is to un- 
derstand the exact role of the impurities. Of course, the 
adjunction of impurities would first decrease the transition 
temperature through the usual dilution effect which, at a 
mean field approximation decreases the average coordina- 
tion number. But in the same time, the number of vor- 
tices (and thus also their disordering consequences) will 
possibly decrease. Also their interactions between each 
other and the interactions between vortices and impuri- 
ties might play some role. This question has been partly 



discussed in refs. [14, 15] and we essentially address here 
the question of the role of impurities at low temperature, 
where the spin-wave approximation should give reliable 
results. 



2 Determination of the phase diagram of the 
diluted model 

The determination of the phase diagram is performed us- 
ing a fit of the order parameter profile inside a finite sys- 
tem to the functional expression predicted by a conve- 
nient conformal mapping, valid at a scale-invariant criti- 
cal point. By extension it is also valid in the whole low- 
temperature phase of the XY model which displays scale- 
invariant algebraic correlation functions. This method has 
been applied with success to the case of the pure XY 
model [16,17,18] and should provide here also reliable re- 
sults. The order parameter vanishes in the bulk of the sys- 
tem at any temperature in the XY model, unless symme- 
try breaking fields are applied along some boundaries dA 
for example. The magnetization profile (cr(r) • (tqa) thus 
obeys a general covariance law under conformal transfor- 
mations. The case of a square geometry with fixed spins 
along its fours edges is particularly easy to implement in 
Monte Carlo simulations. There, the effect of the Schwarz- 
Christoffcl conformal mapping is just to define a rescaled 
distance variable, called k{w) (here and in the following, w 
stands for the complex variable associated to the point r), 
in terms of which one recovers inside the square with fixed 
boundary conditions, a simple power law for the profile: 

{a{w) ■ creA) ^ Hw)]-i^ (5) 

with 



k{'w) = Im 



2Kw' 




2Kw 




X 




r L \ 




L 



Here, cnx, dnx and snx are the Jacobi elliptic functions [19], 
L the linear size of the lattice A and K = K{k) (the com- 
plete elliptic integral of the first kind) and k are constants 
related to the aspect ratio of the system. For more details, 
the reader is referred to rcf. [18]. The main advantage of 
this technique is that one lattice size is in principle suffi- 
cient (provided it is large enough), since the shape effects 
are included in the conformal mapping and the method is 
not much sensitive to finite-size effects. The effect of dis- 
cretization of the lattice is only apparent at the scale of a 
few lattice spacings. One more advantage is the fact that 
all the information encoded in the profile is used, since 
all the points w inside the square arc taken into account 
in the fit. The strategy to obtain the phase diagram is 
summarised in the following: 

i) Perform simulations on a small system (32 x 32) for 
several dilutions (from p = 0.50 to 1.00 every 0.05) and 
several temperatures (from ksT/J = 0.05 to 1.25 every 
0.05). 
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ii) Fit the data to equation (5) (from now on, we use 
the notation m{w) for (cr(u')-crgyi)), compute from a least- 
square fit the sum of deviation squares, x^i per degree of 
freedom (d.o.f.). 

iii) Plot the chi-square per degree of freedom, /d.o.f 
as a function of temperature for all dilutions. It should 
keep a small value (the better the disorder statistics the 
smaller the expected value) in the scale-invariant phase 
and then increase sharply in the paramagnetic phase where 
expression (5) is no longer valid. 

iv) Report on the phase diagram the temperature where 
the /d.o.f . starts to increase sharply, and which is the 
estimate of the critical temperature at a given dilution. 
The results will be refined later when larger lattices will 
be considered in order to investigate the universality class 
at the transition. 

We may also mention here a few comments about the 
technical details: we use a Wolff algorithm [20] and each 
sample is thermalized by 10^ Wolff iterations, and 10'* 
other sweeps are used for computation of the physical 
quantities. The average over disorder realizations is per- 
formed over 10"* samples. One iteration takes of the order 
of 1 /is of CPU time per spin on a standard processor, so 
one simulation needs around 200 x seconds (and here 
we have 23 temperatures times 10 dilutions). For more 
precise estimates to be described later, we also performed 
simulations for sizes 64 x 64 (with 4 • 10'' samples for dis- 
order average, so one simulation takes 8000 x L? seconds) 
and 128 x 128 (with 2 • 10^ samples for disorder average, 
or 4000 X seconds for one simulation). 

In order to illustrate the above mentioned programme, 
we first show in figure 1 the behaviour of the order param- 
eter profile for several temperatures at a dilution p — 0.80 
as a function of the rescaled variable k{w) on a log-log 
scale. It is particularly clear that the data fit quite nicely 
to a power law at low temperatures while this type of fit 
completely breaks down at higher temperatures. 

A 3(i plot of the in the temperature-dilution plane 
is then presented in figure 2. For the sake of clarity, a cut- 
off at 0.5 was introduced in order to avoid too large 
at high temperatures. The low temperature phase with 
quasi-long-range order extends in the whole region in the 
(p, T) plane where the is close to zero, revealing that 
equation (5) nicely fits the numerical data. The phase di- 
agram which can be deduced from these data is shown in 
figure 3 where the data from larger system sizes are also 
reported, as well as previous results [15]. As expected, as 
a result of the infiuence of dilution, the transition temper- 
ature decreases from the pure system value at p = 1 and 
it vanishes at the percolation threshold of site percolation 
on the square lattice. 

A simple standard mean field argument for dilution 
in the low impurity concentration regime [p close to 1) 
gives the beginning of the transition line in the vicin- 
ity of the pure system. The coordination number z on 
the lattice with zero impurity concentration becomes pz 
when only a (small) fraction 1 — p of sites is unoccu- 
pied (this is of course correct only to leading order in 
the very neighbourhood of p = 1, since the problem un- 
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Fig. 1. Order parameter profile for several temperatures 
{ksT/J = 0.1, 0.3; 0.5; 0.7 and 0.9 from top to bottom) at 
a dilution p — 0.80 as a function of the rescaled variable k{w). 
The transition is below ksT/J = 0.5. 



der interest here is site dilution and not bond dilution) 
and thus the transition temperature may be estimated by 
kBTwT{p)/J = pkeTwrip = I)/ J ^ 0.893p. This tran- 
sition curve is shown in dashed line in figure 3 and fits 
correctly only the first point aX p = 0.95. 

In the vicinity of the percolation threshold, one would 
expect that the transition temperature increases with the 
fraction of sites belonging to the percolating cluster. From 
this argument a power-law behaviour would follow, Tc ~ 
(p — Pc)^^'^^7 with the exponent of the order parameter in 
2d percolation (see for example Refs. [21,22]). We cannot 
check this behaviour in the present model where we essen- 
tially focus on smaller dilutions (p closer to 1) and we did 
not perform enough simulations close to the percolation 
threshold. 




k^T/J 



Fig. 2. 3d plot of the x (to get it per d.o.f., one has to divide 
by L^) as a function of temperature and dilution for a system 
of size 32 X 32. 
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Fig. 3. Phase diagram of the 2d dilute XY model. The data 
points correspond to the results obtained after fitting the order 
parameter to the conformal expression. They correspond to the 
temperature where the has a sudden increase. The error bars 
are estimated empirically. Three different sizes are represented 
and the solid line is only a guide for the eyes. The full symbols 
also reported were taken from the work of Leonel et al. [15]. 
Dashed line shows the mean field prediction. 



3 Low temperature behaviour of the diluted 
XY model 

In this section, we address the question of the behaviour 
of the diluted XY model at low temperatures and of the 
universality class along the transition line. More precisely, 
does the spin-wave description provide the correct be- 
haviour at low temperatures in the presence of dilution, 
and does one recover the Kosterlitz-Thouless universal- 
ity class at the transition temperature between the phase 
with quasi-long-range order and the paramagnetic phase? 

First of all, in order to bring reliable arguments, we 
have to produce more precise data at larger sizes to ap- 
proach the thermodynamic limit. We performed simula- 
tions for systems of sizes L = 64 and 128, but only for 
a few temperatures in the vicinity of the transition (es- 
timated from the previous results for the system of size 
L = 32) and for three dilutions, p = 0.70, 0.80 and 0.90. 
The data are shown in figure 4. On the left we plot the 
r] exponent obtained from a least square fit of liim{w) 
vs hiK(w) and on the right, we show the corresponding 
/ d.o.f .. The vertical stripe in each figure is a rough es- 
timate of the transition temperature, where ry takes its 
KT value j and ^ jd.o.f- is at the edge of the plateau re- 
gion (which means that order parameter density profiles 
are still fitted by power laws, i.e. the system is still criti- 
cal). The larger the system size, the smaller the value of 
T where the /d.o.f. increases and the sharper the in- 
crease (note the logarithmic scale on the vertical axis). We 



p 


koT / T 

r^B ^ j 'J 


?7mc 


diluted 
'/SW 


pure 


u.yo 


U.U4 


U.UU (\^) 


U.UUDy 


U.UUD4 




A no 
U.Oo 




O A1 OO 

O.Oioo 


P\ A1 OT 

U.Uiz ( 




n 1 o 
U. iz 




n none 


U.Uiyi 




0.16 




0.0277 




0.92 


0.04 


0.009(2) 


0.0076 


0.0064 




0.08 


0.018(2) 


0.0152 


0.0127 




0.12 


0.027(2) 


0.0227 


0.0191 




0.16 


0.036(2) 


0.0303 


0.0255 


0.88 


0.04 


0.010(2) 


0.0084 


0.0064 




0.08 


0.021(2) 


0.0168 


0.0127 




0.12 


0.032(2) 


0.0251 


0.0191 




0.16 


0.043(2) 


0.0355 


0.0255 



Table 1. Comparison between the Monte Carlo results and 
the spin-wave approximation for the 2d dilute XY model. The 
numerical data were obtained by fitting the density profile in- 
side a square of size L — 128. For comparison, the results of 
the spin- wave approximation of the corresponding pure system 
are also presented. 



also notice that the absolute value of the /d.o.f. in the 
quasi-long-range ordered phase is meaningless, since it is 
strongly dependent on the number of disorder realisations. 
Here the smaller value is obtained for L = 64 where 40000 
samples were produced while only 20000 were realized at 
L = 128. The horizontal scale has been chosen identical 
for all three figures, since it facilitates the comparison and 
makes obvious the role of dilution. 

Let us now discuss the low temperature limit. There, 
a spin- wave calculation (given in the appendix), based on 
the assumption that the spin disorientation remains small 
and allows to expand the cos in the definition of the Hamil- 
tonian, shows that the correlation function exponent is 
simply modified by a "rescaling" of temperature due to 
the presence of the impurities. While the pure model ex- 
ponent is given by 



2ttJ' 



we get in the case of the disordered system 



diluted 

Vsw 



27rJ(l - 2{l-p))' 



(7) 



(8) 



Within the spin-wave approximation, the ?7-exponent still 
varies linearly with temperature at low temperatures, but 
as expected, it increases faster in the disordered model. 
The role of impurities is simply described by an augmented 
effective temperature. We also note that the value of the 
coefiicient, fcB/27rJ(l — 2(1 — p)), is not simply obtained 
from a replacement of the nearest-neighbour coupling J 
by its average pj in the low dilution limit. 

A comparison with numerical results (in table 1) de- 
duced from Monte Carlo simulations shows that the spin- 
wave approximation fits correctly the numerical results at 
very low temperature and low empty sites concentration. 
It is worth noticing that the numerical data were obtained 
on a finite system, and that increasing the system size still 
decreases slightly the resulting exponent (see fig. 4). 
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Fig. 4. Plot of the exponent as a function of the temperature (left) for three different dilutions. The figures on the right give 
the /d.o.f. corresponding to the fit of the order parameter profiles. The vertical stripe materialises the transition temperature 
where the ^ jd.o.j. suddenly increases (logarithmic scale). 
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4 Conclusions 

We have performed an extensive Monte Carlo study of the 

critical behaviour in the quasi-long-range ordered phase of 
the two-dimensional diluted XY model. According to Har- 
ris criterion, one expects weak disorder to be irrelevant at 
the Kosterlitz-Thouless transition where the pure model 
exhibits essential singularities. The numerical results con- 
firm this picture, since at the transition temperature, we 
observe from figure 4 that the ry exponent is compatible 
with the constant value 4 . 
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^ ^ M / X/// XXXXXX-^-^X / //X/ 
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X 1^ / j^^-^ ^ ^ / / / i( / / XX /X 

X |f / ^ X / X X / xx^-^ / / 

// / j //M XX / / // /x/xx>'x / / 

/ //^ ^^i^if/^/f! X X/////X/ / ^ 

////// //////// ////// // ^ ^ 

Fig. 5. Typical configuration for a system of size 32 x 32, 
p = 0.70 at low temperature ksT/ J = 0.05. 

We also obtained the approximate spin-wave solution 
which describes the properties of the diluted model at very 
low temperatures. The presence of non magnetic impuri- 
ties produces an effective increase of the temperature, but 
does not completely suppress the transition before the 
percolation threshold is reached, as it was suspected in 
ref. [15] where a zero-temperature transition was reported 
at p ~ 0.70. The argument reported in ref. [15] was that 
impurities, when located in the vicinity of a pair, produce 
a repulsive interaction between vortices which facilitates 
pair unbinding. This scenario, together with dilution ef- 
fect, enhances the disordering of the system but is not 
sufficient to prevent quasi-long-range ordering at impu- 
rity concentrations above the percolation threshold. As 
an example, a typical configuration at low temperature is 
shown in fig. 5 for a system of size 32 x 32 at p = 0.70. 

In conclusion, one should stress that the behaviour of 
the model in the low temperature critical phase is in full 
agreement with what might be expected from general rel- 
evance arguments. 

Acknowledgement: This work is supported by the French- 
Venezuelan POP program 'Fluides petroliers' and by the French- 
Ukrainian cooperation Dnipro project. The authors gratefully 



acknowledge both programs for their support. Support from 
the CINES under project c20020622309 is also gratefully ac- 
knowledged. B.B. warmly thanks C. Chatelain for his friendly 
technical cooperation and Yu.H. acknowledges useful discus- 
sions with V. Tkachuk. 

Appendix: spin-wave calculation at low tem- 
perature 

In this appendix, we present the result of a spin-wave 
calculation valid in the low temperature limit. First, we 

rewrite the Hamiltonian (4) of the diluted 2d XY model 
for the case of arbitrary short-range ferromagnetic inter- 
action potential J(|r|): 

H=-\Y, - r'l) COs(^r - ^rOCrCr'. (9) 

Here, cos(^r — ^r') stands for a scalar product of two- 
component unit vectors <T(r) directed by angles 6r, Or' (c.f. 

Eq. (2)). Expanding cos{9r — 9r') in (9) for small difference 
in directions of spins one gets the Hamiltonian (9) in the 
spin-wave approximation: 

= H' + Hi. (10) 

The term H' in the Hamiltonian (10) will be further ab- 
sorbed into the energy reference point, whereas to rewrite 
the remaining part Hi we pass to the Fourier-transformed 
angle variables and the interaction potential. Taking the 
periodic boundary conditions we define the Fourier-trans- 
formed quantities as: 

■^w = ^Ee*'"^('^)' Kk) = Ee"*""-^w- (12) 

k r 

In (11), (12) J2r spans A'' sites of the lattice whereas 
is within the first Brillouin zone. Written in the Fourier- 
transformed variables the Hamiltonian Hi reads: 

^1= I Yl S(q + ki+k2)5(-q)i/(q)^k,^k. 

q,ki,k2 

E ^(q + ki)5(k2-q)Kq)^k,ek„ (13) 

q,ki,k2 

with 

^(q) = ^E^^''""^- (14) 

r 

Let us single out in (13) the Hamiltonian of the pure 
2d XY model. To this end let us introduce variables p(q) 
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which may serve to show deviation of S'(q) from the Kro- 
neker's delta <5(q): 



with 



S(q) = 5(q) - p(q), 



(15) 



(16) 



Denoting a configurational average with the distribu- 
tion function (3) by (....) we get for the two first moments 
of the random variable Cr- = X]cr=o i ^(cr)cr — P and 
CrCr' = for r ^ r', or p for r = r'. Here, we remember 
that p is the concentration of the occupied sites cv = 1 and 
that V[cr] = Ylr ^(cr)- This leads to the following values 
for the first two moments of p(q) (note, that due to (9) r 

r'):'^= (1 -p) 6{q) and p(q)p(q') = {l-pf (5(q)5(q'). 
Considering the concentration of empty sites (1 —p) to be 
small wc take it as an expansion parameter and further 
we keep only linear in p(q) contributions. For the Hamil- 
tonian (9) we get: 



ki,k2 

-Z/(k2)] /9(ki + k2)^k,^k. + O(P'). (17) 



Where Hq stands for the Hamiltonian of the undiluted 
system: 



(18) 



For the nearest neighbour interaction and square lattice 
with a lattice constant a = |/x| we have J(r) = J for 
|r| = a and for |r| ^ a, 



z/(k) = 2J(cos(A;xa)-|-cos(fc3;a)) ~ 4:J-J\'k\'^a'^+. 
Subsequently, the Hamiltonian (17) reads: 

k 

+ ^ Jo2(ki.k2)p(ki+k2)0k,ek., 

ki,k2 



(19) 



(20) 



with a scalar product ki.k2. 

For a given configuration of disorder, the configuratio- 
nally-dependent partition function Zgonf is defined by: 



-'conf 



= Sp^e 



-I3H 



(21) 



where Spg means integration over spin degrees of freedom 
on each site (and /? = (fc^T)"^): 



sp.(...)=n£^(..)- 



(22) 



For a given configuration of occupied and empty sites, let 
us define a thermodynamic averaging by: 



((•••)> = 



^conf 



Sp.e-^^(...). 



(23) 



Considering the quenched dilution (i.e. the case when mag- 
netic and non-magnetic sites arc fixed on their places) one 
gets the observables by averaging thermodynamically av- 
eraged quantities with respect to different configurations 
of disorder [23] . Now, the pair correlation function is de- 
fined as 



G2(|r2 - ri|) = (Cr^Cri C0s(^r2 " ^n)) 



^ -J—Sp,e-0HCr,Cr, (l - hOr, - er,r). (24) 

^conf \ Z J 

As far as we keep only linear in p(q) terms, we may de- 
couple configurational averaging in (24) and write the pair 
correlation function as: 



G2(|r2-ri|) 



Spge ^^Cr2Cri(^l 



Sp,e-/3^^ 



(25) 

For the configurationally averaged value of the Hamilto- 
nian R one gets: 

^5Zj|k|2a2M_k+ J'«'ki.k2p(ki+k2)eki^k2 



ki,k2 



= i^Jjkpa^^k^.k 



+ ^ Ja2ki.k2^(ki + k2)(l - p)eki^k2 



ki,k2 



l[i-2(i-p)]^j|k|VM-k 



[l-2(l-p)]/fo- 



(26) 



Then the formula for the pair correlation function (25) 
reads: 

G2(|r2 - ri|) = const x 

xSp.e-'^'^o^i _ (Q^^ _ ^,j72)/Sp,e-'5'«°, (27) 



with 



const = 1 — 2(1 — p), 
/3' = [l-2(l-p)]/3. 



(28) 
(29) 



Formula (27) for the pair correlation function of the 
diluted 2d XY model has the form of the spin-wave— 
approximated pair correlation function of the pure system. 
It differs only by a constant (28) which renormalizes also 
the temperature f3 (29). Taking that in the spin- wave ap- 
proximation the pair correlation function critical exponent 
for the pure system equals [7]: 



(30) 



one gets the exponent for the diluted system substituting 
in (30) temperature /? by /?' (29): 



^diluted 

Vsw 



1 



27rJ/3(l-2(l-p))' 



(31) 
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